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Summary 


Two  computational  methods  are  considered  for  a class  of  non- 
uniform,  continuous,  multiple  time-delay  linear  systems,  which  have 
been  developed  in  the  modeling  of  lossless  layered  media.  In  the 
first  approach,  we  discretize  the  time  axis  and  insert  states  of 
intermediate  delays,  to  arrive  at  a set  of  standard  finite-difference 
equations.  For  our  particular  system,  matrix  multiplications  can  be 
reduced  to  simple  scalar  multiplications.  In  the  second  approach, 
we  define  mapping  rules  for  the  transformation  of  states  at  an  inter- 
faces, and  keep  a state  reference  table  for  look-up  and  branching. 

The  procedure  is  similar  to  ray-tracing.  Several  experimental  re- 
sults are  presented  to  show  the  trade-off  between  storage  requirement 
and  CPU  time -spent  for  the  two  methods. 


1.  Introduction 


In  the  development  of  time-domain  state  space  models  for  loss- 
less layered  media,  Mendel  et.  al.  [1]  have  arrived  at  a class  of 
non-uniform,  multiple,  continuous,  time-delay  linear  equations. 

These  equations  are  termed  Causal  Functional  Equations.  In  this 
paper  we  will  discuss  two  computational  methods  to  obtain  a solution 
for  this  class  of  equations. 

The  immediate  application  of  this  solution  is  to  generate 
synthetic  seismograms  for  geophysical  study  (e.  g. , in  testing  decon- 
volution algorithms  and  inversion  techniques).  As  mentioned  in  the 
paper  by  Mendel  et.  al.  [1],  this  class  of  equations  can  also  be  used 
to  model  lossless  transmission  lines  and  certain  acoustic,  optical 
and  electromagnetic  systems. 

In  the  next  section  we  will  give  a formulation  of  this  class  of 
equations  and  discuss  some  of  its  characteristics.  In  Section  3 we 
discuss  the  method  of  discretization  of  the  time-axis  and  the  proce- 
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dure  to  obtain  a standard  unit- 
delay  finite-difference  equations. 

A ray  tracing  technique  using 
mapping  rules  for  states  is  dis- 
cussed in  Section  4.  In  Section 
5 we  make  a comparison  between 
the  two  methods  in  terms  of  sim- 
plicity, storage  requirement  and  CPU 
time-spent.  Some  experimental  re- 
sults are  presented.  We  conclude 
with  a summary  of  results  and  a sam- 
ple plot  of  unit  pulse  response  for 
layered  earth  media. 

2.  Causal  Functional  Equations 

A system  of  K layered 
media  is  depicted  in  Figure  1. 

For  each  layer,  we  define  two 
states:  uj(t)  and  dj(t),  which 

are  the  upgoing  ana  downgoing 
states  at  the  (j-l)th  and  jth 

interfaces  respectively.  The  reflection  and  transmission  characteris- 
tics at  an  interface  are  modeled  by  the  reflection  coefficients,  r-  . 
Within  the  layer,  the  state  only  suffer  a time  delay,  tj  . A more 
detailed  description  can  be  found  in  Mendel  et.  al.  [1], 


~Tj7rrTjrj-T — ,K-2’'W 

Zffpnyyyf ? 


'.nttrfoct  - * l r_ ) 


'N  ^ ' '"'N  - V/'v 

\ \ 

5 f f / ) • / 1 i / 

) / / j ;/ )>  f n 


Figure  1.  System  of  K layered 
media. 


The  Causal  Functional  Equations  are: 


where 


d(t  + t) 

= Aj  d(t)  + A2  u(t)  + £ m(t) 

( 1 ) 

U(t  + T) 

= A3  d(t)  + A4  u(t) 

( 2 ) 

z(t) 

= h u(t)  + rQ  m(t) 

( 3 ) 

d(t+jr)  = coMd^t  + Tj),  d2(t+ t2),  . . . ,dK(t  + tk)) 

U(t+T)  = COl(ui(t  + T1),U2(t  + T2),  . ...U^t+1^)) 
d(t)  = col(dj(t),  d 2(t) , ....  d^(t)) 
u(t)  = coltu^t),  u2(t),  . . . , uK(t)) 


( 4 ) 


and 


Aj  = sub-diag  (1+^,1+^ l + rK1) 


( 5 ) 


( 6 ) 
( 7 ) 


A3  = diag  (r  j , r ^ rR  ) 

A4  = super  -diag  (l-r^.l-r^ ,1_rK-l^ 

£ = col  (1  + r , 0,  0,  . . . , 0) 
h = col  (1  - rQ  , 0,  0, . . . , 0) 

The  notation  is: 


( 8 ) 

( 9 ) 

( 10  ) 


with  the  non-zero  entries  given  in  the  above  equations.  Matrices  A^  , 
A2,A3'A4  are  K*K>  Sl>  h are  KX  1.  Physically,  Eqs.  (l)-(3)  des- 
cribe the  interactions  and  transformations  of  the  downgoing  and  up- 
going  states  in  a layered  media.  These  equations  are  continuous- 
time equations  with  non-uniform,  multiple  time-delays.  They  are 
linear  in  the  states;  hence,  we  may  decompose  them  in  suitable  form 
and  apply  the  principle  of  superposition  to  obtain  a solution.  The 
output  response,  y(t),  can  be  obtained  by  the  convolution  of  the  input 
wave  form,  m(t),  and  the  impulse  response  of  the  system,  which  is 
itself  a sequence  of  non-uniformly  spaced  impulses  along  the  time- 
axis.  In  general,  those  impulses  will  become  more  densely  populated 
on  the  time-axis  as  time  increases. 


3.  Finite-Difference  Technique 


The  system  of  equations  in  Section  2 can  be  converted  into 
standard  finite-difference  equations  (FDE)  by  discretization  of  time. 
Since  the  delay  terms  are  non-uniform  and  not  multiple  of  one  an- 
other, we  may  have  to  choose  a very  small  time  interval,  A,  for 
discretization.  Below,  we  show  the  procedure  to  arrive  at  the  de- 
sired FDE; 

First,  let  us  find  A,  the  greatest  common  factor  of  the  delay 
terms,  such  that  Tj  = k^A , i=l,2,  ...,K.  The  kj^  are  positive  integers. 
Let  t = kA  and  define  x(k)=x(kA).  We  make  the  following  transformations 
in  Equ.  (1)  and  (2); 

d(t+Tj  = d(k  + k^) 

u(t  + T)  = u(k+k.) 

Next,  for  each  state  with  ki>l,  we  insert  states  of  intermediate 
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delays  to  obtain  a unit-delay  system.  Specifically,  we  make  the 
following  changes; 

d.(k  + k.)  -*  col  (d.fk  + k.bd.fk+k.-  1) ,d.(k+l)),  i=  1,2 K. 

u^k  + k.)  -»  col  (u.(k+  1),  u.(k  + 2),  , u^k  +k.)),  i=  1,  2,  . . . , K. 

The  above  ordering  of  dj(k)  and  uj(k)  preserves  the  structure  of  the 
matrices  Aj.A2.A3  and  A4  . After  staking  up  the  states  into  vectors 
d(k+l)  and  ja(k+l),  we  arrive  at  the  following  equations: 

. d(k  + 1 ) = A j d(k)  + A ^ u(k)  + £1  m(k)  ( 11  ) 

u(k+  1)  = A*  d(k)  + A*  u(k)  ( 12  ) 

£(k)  = h1  u(k)  + rQ  m(k)  ( 13  ) 

where 

d(k+l)  k col  (djfk  + kj). d^k+l)  j d2(k  + k2),  ...,d2(k+l)  | ...  | 

dK(k  + kK)"“,dK(k+ 1)J  ( 14  ) 

u(k+l)  = col  (u1(k+l),...fu1(k+k1)  j u2(k+l), u2(k  + k2)  | ...  ! 

uK(k+  1 ),...,  u^(k  + kK>)  (15) 

and 

Aj  = sub-diag  (1,  1 1 j (1  + ri),  1,  ....  1 j . ..  j (1  + rK_ j )»  1 1) 

( 16  ) 

A2  = -diag  (rQ  , 0, . . . , 0 j r^  0,  ....  0 j . . . j r^,  0,  ....  0)  ( 17  ) 

a!  = diag  (0 ,0,  r.  ! 0 0,  r.  | ...  J 0 ,0,  r ) ( 18  ) 

J 1 I £ 1 1 K. 

A4  = super -diag  (1,  . . . , 1,  (1  - r ) J . . . { 1. ....  1,  (1  - r ) ! 1 1) 

* 1 i I k.  - 1 1 

( 19  ) 

£ = col  (1  + rQ  , 0,  ....  0 j 0, ....  0 j ...  j 0, ....  0)  ( 20  ) 

hl  = col  (1  - rQ,0, . . .,  0 j 0, ...  ,0  j . . . j 0,  ....  0)  (21) 

Let  L=  £ k.  , then,  d(k  ),u(k)  , g^h1  are  Lx  1,  a!  , A*  , A*  , A*  are 
j=  j 1 — — — 1234 

Lx  L. 


To  obtain  the  system's 
unit  pulse  response  we  set 
m(0)=l,  and  m(k)  = 0,  ki  1. 
Our  computational  algorithm 
for  observation  time  up  to 
TMAX  is: 

(1)  Initialization:  k=0 
d(l)  = s!  ) 

u(l)  = 0 j ( 22  ) 

y(°)  = ^ 


jd,(K*2) 
_ Jd,(KM) 
jd|IK) 


V(K) 


’ WT 


u,(K*  2) 


-InteHoce-i  I r ) 


— (r 1 0) 
& 

— Ir  *0) 


-Inierlace-i ( r.) 


Figure  2.  Illustration  to  show  in- 
sertion of  new  state  vari- 
able for  the  ith  layer 
where  k.  = 3A  . 


(2)  Recursive  Step:  k=l,2, 


TMAX 


d(k+  1)  = aJ  d(k)  + A*  u(k) 
u(k  + 1)  = A*  d(k)  + A1a  u(k ) 
y(k)  = h1  u(k) 

Repeat. 


( 23  ) 


Since  each  row  of  matrices  Aj.A2.A3  and  A4  has  only  one 
non-zero  element,  the  multiplication  indicated  in  Eq.  (23)  is  just  a 
collection  of  scalar  multiplications.  The  non-zero  entries  of  the 


four  Aj  matrices  are  related  and  can  be  obtained  from  a single 
array.  A,  with  elements  (r0:  A^  ).  New  values  of  u(k)  and  d(k)  can 
be  stored  back  into  the  state  vectors;  thus,  we  only  need  three  arrays 
for  A,  u(k),  and  d(k). 


For  our  particular  system,  the  insertion  of  states  corresponds 
to  the  dividing  of  a layer  into  equal  sub-layers  by  inserting  inter- 
faces whose  reflection  coefficients  are  zero.  This  phenomenon  can 
be  observed  from  the  A1  matrices  in  Eqs.  ( 1 6) - ( 1 9).  The  particular 
ordering  of  the  states  u(k)  and  d_(k)  can  be  understood  from  the  way 
the  states  are  defined  and  is  illustrated  for  the  case  kj=  3a  in  Figure  2. 

4.  Ray  Tracing  Technique 

In  this  section,  we  utilize  the  properties  of  linearity  and  the 
causal  nature  of  our  equations  to  develop  a ray  tracing  technique. 
Instead  of  focusing  on  the  interaction  of  states,  we  decompose  Eqs. 

( 1 ) - ( 3 ) and  define  the  following  mapping  rules  to  track  how  a state 
propagates  at  an  interface: 


( 24  ) 


d (t) 


u (t) 


m(t)  -* 

y(t) 

j = 0 

ri  <Vt| 

j = 1 9 • • 

• » K 

(l+r.)d.(t)  -* 

J J 

j - 0,  . . 

• * K- 

(i  - r0}  ui(t>  - 

y(t) 

j = 1 

(1  - r ) u (t)  -* 

uj-i<,+Tj-i| 

j = 2,  . . 

• * K 

-ri-i  “j*"  ^ 

d.(t  + T.| 

j = 1 , . . 

• » K 

( 25  ) 


We  illustrate  the  interpretation  of  these  mapping  rules  for 
d.(t),  in  Eq.  (24)  in  Figure  3. 
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Figure  3.  The  transformation  of  d:(t)  into  neighboring  states  at 

(a)  the  surface,  (b)  jth  interface,  and  (c)  bottom  interface. 

In  order  to  use  Eqs.  (24)  and  (25),  we  need  to  have  a state 
reference  table.  Each  element  of  the  states  in  the  table,  called  an 
event,  is  characterized  by  its  time  and  amplitude.  A subroutine, 
ISMAL,  picks  up  the  event  of  the  next  smallest  time  and  processes 
it  according  to  the  mapping  rules,  which  are  coded  into  the  UP  and 
DOWN  subroutines.  The  new  events  generated  are  stored  back  into 
the  table.  As  each  event  branches  into  two  new  events,  the  table 
grows  geometrically.  To  restrict  the  growth  of  the  table,  we  use 
two  tolerance  parameters,  AMIN  and  6T,  and  collapse  events  that 
occur  at  the  same  time  point  before  we  store  them  back  into  the 
table.  AMIN  controls  the  amplitude  of  the  computed  event.  If  that 
amplitude  is  less  than  AMIN,  it  is  set  to  zero.  6T  controls  the  time 
separation  of  two  events  below  which  they  are  considered  to  occur  at 
the  same  time.  The  output  is  observed  up  to  a pr e-determined  time, 
TMAX. 


Comparison  of  the  Two  Methods 


If  we  expand  the  matrix  multiplications  in  the  Section  3 method, 
we  obtain  the  Section  4 mapping  rules,  from  step  k to  k+1;  hence, 
the  two  methods  are  related.  In  the  discrete  method  we  keep  the 
history  in  the  state  vectors,  and  in  the  ray-tracing  technique  we  use 
a state  reference  table. 

The  discrete  method  is  simple,  easy  to  implement  and  recur- 
sive; however,  due  to  the  non-uniform  nature  of  the  delay-terms,  we 
may  have  to  use  a A small  enough  such  that  it  is  a submultiple  of 
all  the  delay  terms.  The  storage  requirement  is  usually  not  exces- 
sive. However,  we  may  be  computing  zero  for  a lot  of  points  on 
the  time  axis  where  no  events  occur.  Thus  the  CPU  time  may  be 
larger  than  for  the  ray-tracing  technique. 

The  ray-tracing  technique  requires  a slightly  more  complex 
program  to  implement  it.  Its  major  disadvantage  is  the  huge  storage 
requirement  for  the  state  reference  table.  Since  this  technique  com- 
putes only  the  actual  event  points,  it  can  be  very  efficient  and  fast 
for  short  observation  time. 

We  present  below  three  experimental  results  in  testing  the  two 
methods,  using  a PDP-10  KI  interactive  system: 

Experiment  1:  3 layers,  TMAX=5.00  minutes,  A = 0.001  minute, 

AMIN  = 10*^0,  6T  = 10“^  minutes. 

t = (0.3,  0.001,  0.5)  minutes 
r = (0.8,  -.3,  .3,  .5) 

This  experiment  illustrates  the  potential  problem  of  computing  arti- 
ficial zero  points  for  FDE.  The  first  non-zero  data  occurs  at 
k=  600  1 

Experiment  2:  3 layers,  TMAX=10  minutes,  A = 0.01  minute, 

AMIN  = 10'  1 , 6T=10"5  minutes. 

T = (0.03,  0.01,  0.  05)  minutes 
r = (0.8,  -.3,  .3,  .5) 

This  experiment  illustrates  storage  problem  for  ray  tracing  technique. 

Experiment  3:  19  layers,  TMAX  = 2.  20  minutes,  A = 0.001, 

AMIN  = 1 0 ' 0 , 6T  = 10“^  minutes. 


T 


(.016, 

.050, 

.004, 

.023, 

.022, 

.015,  . 

042,  . 028,  .006, 

. 038, 

. 003, 

. 006, 

. 007, 

. 072, 

. 005,  . 

030,  .027,  .038, 

.014) 

minutes 

(.99,  , 

.028,  - 

.061, 

. 082, 

. 034, 

-.  068, 

-.  016,  .168,  -.008, 

00 

o 

.058, 

. 114, 

-. 057, 

. 026. 

-.112, 

-.220,  .076,  .156, 

. 039, 

-. 229) 

This  experiment  illustrates  the  problems  of  CPU  and  storage  for  a 
large  number  of  layers  case. 


Experiments 

FDE 

Ray  Tracing 

(1)  CPU 

4 min.  46.  05  sec. 

4.70  sec. 

Storage 

7400 

8726 

Number  of  Output 

5000 

436 

(2)  CPU 

2.  35  sec. 

2.  85  sec. 

Storage 

1027 

3434 

Number  of  Output 

1000 

163 

(3)  CPU 

1 min.  20.  80  sec. 

1 min.  4.  56  sec. 

Storage 

3538 

73902 

Number  of  Output 

2200 

976 

6.  Conclusions 


In  this  paper  we  have  presented  two  computational  methods  for 
a class  of  causal  functional  equations.  They  are  the  discrete  method 
of  finite-difference  equations  and  the  continuous  method  of  ray-tracing 
The  FDE  has  low  storage  requirement  but  may  run  into  high  CPU 
time-spent  if  the  sampling  interval.  A,  is  too  small.  The  ray 
tracing  technique  has  a large  storage  requirement  but  can  be  very 
efficient  for  short  observation  time  and  small  number  of  layers  cases 
The  actual  choice  of  a particular  method  will  depend  on  the  input 
data  set.  A sample  plot  of  the  unit  impulse  response  for  Experi- 
ments 1 and  2 is  illustrated  in  Figure  4. 
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Figure  4.  Unit  Impulse  Response  for  Experiments  I & 2.  Observe  the 
rapid  decay  in  amplitudes  of  the  impulses. 
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